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Abstract — This paper addresses the problem of expressing a 
signal as a sum of frequency components (sinusoids) wherein 
each sinusoid may exhibit abrupt changes in its amplitude 
and/or phase. The Fourier transform of a narrow-band signal, 
with a discontinuous amplitude and/or phase function, exhibits 
spectral and temporal spreading. The proposed method aims 
to avoid such spreading by explicitly modeling the signal of 
interest as a sum of sinusoids with time-varying amplitudes. 
So as to accommodate abrupt changes, it is further assumed 
that the amplitude/phase functions are approximately piecewise 
constant (i.e., their time-derivatives are sparse). The proposed 
method is based on a convex variational (optimization) approach 
wherein the total variation (TV) of the amplitude functions are 
regularized subject to a perfect (or approximate) reconstruction 
constraint. A computationally efficient algorithm is derived based 
on convex optimization techniques. The proposed technique 
can be used to perform band-pass filtering that is relatively 
insensitive to narrow-band amplitude/phase jumps present in 
data, which normally pose a challenge (due to transients, leakage, 
etc.). The method is illustrated using both synthetic signals and 
human EEG data for the purpose of band-pass filtering and the 
estimation of phase synchrony indexes. 

Index Terms — sparse signal representation, total variation, 
discrete Fourier transform (DFT), instantaneous frequency, phase 
locking value, phase synchrony. 

I. Introduction 

SEVERAL methods in time series analysis aim to quantify 
the phase behavior of one or more signals (e.g., studies 
of phase synchrony and coherence among EEG channels 
[5], [6], [28], [33], [64], [65]). These methods are most 
meaningfully applied to signals that consist primarily of a 
single narrow-band component. However, in practice, available 
data often does not have a frequency spectrum localized to 
a narrow band, in which case the data is usually band-pass 
filtered to isolate the frequency band of interest (e.g., [33], 
[64]). Yet, the process of band-pass filtering has the effect of 
spreading abrupt changes in phase and ampUtude across both 
time and frequency. Abrupt changes present in the underlying 
component will be reduced. This limitation of linear time- 
invariant (LTI) filtering is well recognized. In linear signal 
analysis, time and frequency resolution can be traded-off with 
one another (c.f. wavelet transforms); however, they can not 
both be improved arbitrarily. Circumventing this limitation 
requires some form of nonlinear signal analysis [14J. The 
nonlinear analysis method developed in this paper is motivated 
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by the problem of extracting (isolating) a narrow-band signal 
from a generic signal while preserving abrupt phase-shifts and 
amplitude step-changes due for example to phase-resetting 
[56]. 

This paper addresses the problem of expressing a signal 
as a sum of frequency components (sinusoids) wherein each 
sinusoid may exhibit abrupt changes in its amplitude and/or 
phase. The discrete Fourier transform (DFT) and Fourier series 
each give a representation of a generic signal as a sum of 
sinusoids, wherein the amplitude and phase of each sinusoid 
is a constant value (not time- varying). In contrast, in the 
proposed method, the amplitude and phase of each sinusoid is 
a time-varying function. So as to effectively represent abrupt 
changes, it is assumed that the amplitude and phase of each 
sinusoid is approximately constant; i.e., the time-derivative 
of the amplitude and phase of each sinusoid is sparse. It is 
further assumed that the signal under analysis admits a sparse 
frequency representation; i.e., the signal consists of relatively 
few frequency components (albeit with time-varying amphtudes 
and phases). 

In the proposed sparse frequency analysis (SFA) approach, 
the amphtude and phase of each sinusoid is allowed to vary 
so as (0 to match the behavior of data wherein it is known or 
expected that abrupt changes in amplitude and phase may be 
present (e.g., signals with phase-reset phenomena), and {ii) to 
obtain a more sparse signal representation relative to the DFT 
or overs ampled DFT, thereby improving frequency and phase 
resolution. 

The proposed method has a parameter by which one can 
tune the behavior of the decomposition. At one extreme, the 
method coincides with the DFT. Hence, the method can be 
interpreted as a generalization of the DFT. 

In order to achieve the above-described signal decomposition, 
we formulate it as the solution to a convex sparsity-regularized 
linear inverse problem. Specifically, the total variation (TV) of 
the amplitude of each sinusoid is regularized. Two forms of the 
problem are described. In the first form, perfect reconstruction 
is enforced; in the second form, the energy of the residual 
is minimized. The second form is suitable when the given 
data is contaminated by additive noise. The two forms are 
analogous to basis pursuit (BP) and basis pursuit denoising 
(BPD) respectively [14]. 

To solve the formulated optimization problems, iterative 
algorithms are derived using convex optimization techniques: 
variable splitting, the alternating direction method of multipliers 
(ADMM) [3], [7], [15], and majorization-minimization (MM) 
[21]. The developed algorithms use total variation denoising 
(TVD) [12], [54] as a sub-step, which is solved exactly 
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using the recently developed algorithm by Condat [16]. The 
resulting algorithm is 'matrix-free' in that it does not require 
solving systems of linear equations nor accessing individual 
rows/columns of matrices. 

The proposed signal decomposition can be used to perform 
mildly nonlinear band-pass filtering that is better able to track 
abrupt amplitude/phase jumps than an LTI bandpass filter. Such 
band-pass filtering is also relatively insensitive to narrow-band 
amplitude/phase jumps present outside the frequency band of 
interest (interference), which normally pose a challenge (due to 
transients, leakage, etc.). The method is illustrated using both 
synthetic signals and human EEG data, for the purpose of band- 
pass filtering and the estimation of phase synchrony indexes. 
In the examples, the proposed method is compared with the 
DFT and with band-pass filtering. The examples demonstrate 
the improved ability to represent and detect abrupt phase shifts 
and amplitude changes when they are present in the data. 

A. Related work 

The short-time Fourier transform (STFT) can already be 
used to some extent to estimate and track the instantaneous 
frequency and amphtude of narrow-band components of a 
generic signal [4], [23], [34]; however, the STFT is computed 
as a linear transform and is defined in terms of pre-defined basis 
functions, whereas the proposed representation is computed 
as the solution to a non-quadratic optimization problem. As 
discussed in [48], the time-frequency resolution of the STFT is 
constrained by the utihzed window. Hence, although the STFT 
can be computed with far less computation, it does not have 
the properties of the proposed approach. On the other hand, the 
effectiveness and suitabihty of the proposed approach depends 
on the vahdity of the sparsity assumption. 

We also note that the sparse STFT (i.e. basis pursuit with 
the STFT) and other sparse time-frequency distributions can 
overcome some of the time-frequency Umitations of the lineeir 
STFT [22]. However, the aim therein is an accurate high 
resolution time-frequency distribution, whereas the aim of this 
paper is to produce a representation in the time-domain as a 
sum of sinusoids with abrupt amplitude/phase shifts. 

Sinusoidal modeling also seeks to represent a signal as 
a sum of sinusoids, each with time- varying amplitude and 
phase/frequency [35], [37], [43], [44]. However, sinusoidal 
models usually assume that these are slowly varying functions. 
Moreover, these functions are often parameterized and the 
parameters estimated through nonlinear least squares or by 
short-time Fourier analysis. The goal of this paper is somewhat 
similar to sinusoidal modeling, yet the model and optimization 
approach are quite distinct. The method presented here is non- 
parametric and can be understood as a generalization of the 
DFT. 

Empirical mode decomposition (EMD) [30], a nonlinear data- 
adaptive approach that decomposes a non-stationary signal into 
oscillatory components, also seeks to overcome limitations of 
linear short-time Fourier analysis [1], [51], [52]. In EMD, the 
oscillatory components, or 'intrinsic mode functions' (IMFs), 
are not restricted to occupy distinct bands. Hence, EMD 
provides a form of time-frequency analysis. In contrast, the 



sparse frequency analysis (SFA) approach presented here is 
more restrictive in that it assumes the narrow band components 
occupy non-overlapping frequency bands. EMD, however, is 
susceptible to a 'mode mixing' issue, wherein the instantaneous 
frequency of an IMF does not accurately track the frequency 
of the underlying non- stationary component. This is the case, 
in particular, when the instantaneous amplitude/phase of a 
component exhibits an abrupt shift. Several optimization-based 
variations and extensions of EMD have been formulated that 
aim to provide a more robust and flexible form of EMD 
[27], [39], [40], [45], [49]. In these methods, as in EMD, 
one narrow-band component (IMF) is extracted at a time, the 
frequency support of each component is not constrained a priori, 
and the instantaneous amplitude/phase of each component 
is modeled as slowly varying. In contrast, the SFA method 
developed below, optimizes all components jointly, uses a 
fixed vmiform discretization of frequency, and models the 
instantaneous amplitude/phase as possessing discontinuities. 
Hence, SFA and EMD-like methods have somewhat different 
aims and are most suitable for different types of signals. 

Other related works include the synchrosqueezed wavelet 
transform [17], the iterated Hilbert transform (IHT) [24], 
and more generally, algorithms for multicomponent AM-FM 
signal representation (see [25]). These works, again, model the 
instantaneous amplitudes/phase fimctions as smooth. 



II. Preliminaries 

A. Notation 

In this paper, vectors and matrices are represented in bold 
(e.g. X and D). Scalars are represented in regular font, e.g., A 
and a. 

The VV-point sequence x{n), n G Z^r = {0, . . . , VV — 1} is 
denoted as the column vector 

^=[x{0), x{l), x{N-l)]\ (1) 

The A'^-point inverse DFT (IDFT) is given by 

fe=0 

The IDFT can alternately be expressed in terms of sines and 
cosines as 

x{n) = ^ {akCOs(^^^nj + bksm(^^^njy (3) 

fc=0 

We use k to denote the discrete frequency index. The 
variables c{n, k) and s(n, k) will denote c(n, k) = cos{2tt fku) 
and s(n, k) = sin(27r/fcn) for n G Zjv and k G Zk = 
{0,...,K~l}. 

When a s M.^^^ is an array, we will denote the A;-th 
colunm by afe, i.e., 

a= [ao, ai, a^-i]- (4) 

The ii norm of a vector v is denoted ||v||i = The 
energy of a vector v is denoted ||v||2 = 
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B. Total Variation Denoising 

The total variation (TV) of tiie discrete A/'-point signal x is 
defined as 

N-l 

TV(x)= ^ \x{n)-x{n-l)\. (5) 

n=l 

The TV of a signal is a measure of its temporal fluctuation. It 
can be expressed as 



TV(x) = ||Dx||i 



where D is the matrix 



D = 



-1 1 



-1 1 



-1 1 



(6) 



(7) 



of size (N — 1) X N. That is, Dx denotes the first-order 
difference of signal x. 

Total Variation Denoising (TVD) [54] is a nonlinear filtering 
method defined by the convex optimization problem, 



tvd(y, A) := argniin||y - x||2 + A||Dx||i 



(8) 



where y is the noisy signal and A > is a regularization 
parameter. We denote the output of the TVD denoising 
problem by tvd(y, A) as in (8). Specifying a large value for 
the regularization parameter A increases the tendency of the 
denoised signal to be piece-wise constant. 

Efficient algorithms for ID and multidimensional TV de- 
noising have been developed [11], [13], [42], [53], [61] and 
applied to several linear inverse problems arising in signal and 
image procesing [19], [20], [41], [55], [62]. Recently, a fast 
and exact algorithm has been developed [16]. 

III. Sparse Frequency Analysis 

As noted in the Introduction, we model the A^-point signal 
of interest, x e M^, as a sum of sinusoids, each with a 
time-varying amplitude and phase. Equivalently, we can write 
the signal as a sum of sines and cosines, each with a time- 
varying amplitude and zero phase. Hence, we have the signal 
representation: 

K-l 

= (a/cH cos(27r/fen) + bk{n) sin(27r/fen)^ (9) 

where //,, = k/K and ak{n), hk{n) € M, fc G Z^, n £ Z^. 

The expansion (9) resembles the real form of the inverse 
DFT (3). However, in (9) the amplitudes Ofc and bk are time- 
varying. As a consequence, there are 2NK coefficients in (9), 
in contrast to 2N coefficients in (3). 

In addition, note that the IDFT is based on N frequencies, 
where N is the signal length. In contrast, the number of 
frequencies K in (9) can be less than the signal length N. 
Because there are 2KN independent coefficients in (9), the 
signal x can be perfectly represented with fewer than N 
frequencies (K < N). 

In fact, with if = 1, we can take ao{n) = x{n) and bo{n) 
I) and satisfy (9) with equality. [With = 1, the only term 



in the sum is ao(n) cos(27r/on) with /o = 0. ] However, 
this solution is uninteresting - the aim being to decompose 
x{n) into a set of sinusoids (and cosines) with approximately 
piecewise constant ampUtudes. 

With K > 1, there are infinitely many solutions to (9). 
In order to obtain amplitude functions afe(n) and &ft(n) that 
are approximately piecewise constant, we regularize the total 
variation of G and G M.^ for each k G Zk, where 
Bfe denotes the A^-point amphtude sequence of the fc-th cosine 
component, i.e., afc = ((ifc(n))„gZN- Similarly, bfc denotes 
the iV-point amplitude sequence of the fc-th sine component. 
Regularizing TV(afe) and TV(bfe) promotes sparsity of the 
derivatives (first-order difference) of afc and b^, as is well 
known. Note that TY{ak) can be written as ||Dafc||i. This 
notation will be used below. 

It is not sufficient to regularize TV(afc) and TV(bfc) 
only. Note that if afe(n) and bk{n) are not time-varying, 
i.e., afc(n) = Uk and bk{n) = bk, then TV(afc) = and 
TV(bfc) = 0. Moreover, a non-time-varying solution of this 
form in fact exists and is given by the DFT coefficients, c.f. 
(3). This is true, provided K ^ N. That is, based on the DFT, 
a set of constant-amplitude sinusoids can always be found so 
as to perfectly represent the signal x{n). The disadvantage 
of such a solution, as noted in the Introduction, is that if 
x{n) contains a narrow-band component with amplitude/phase 
discontinuities (in addition to other components), then many 
frequency components are need for the representation of the 
narrow-band component; i.e., its energy is spectrally spread. 
(In tum, due to each sinusoid having a constant amplitude, the 
energy is temporally spread as well.) In this case, the narrow- 
band component can not be well isolated (extracted) from the 
signal x(n) for the purpose of further processing or analysis 
(e.g., phase synchronization index estimation). 

Therefore, in addition to regularizing only the total variation 
of the amplitudes, we also regularize the frequency spectrum 
corresponding to the representation (9). We define the frequency 
spectrum, Zk, as 



n 

or equivalentiy 



|6fe(n)|2, keZK 



afc 



li + l|bfc||i, 



k e 



(10) 



(11) 



The frequency-indexed sequence Zk measures the distribution 
of amplitude energy as a function of frequency. Note that, in 
case the amphtudes are not time- varying, i.e., ak{n) = cik and 
bk{n.) = hk, then \zk\ represents the modulus of the fc-th DFT 
coefficient. We define the vector z G as z = {zk)keZK- In 
order to induce z to be sparse, we regularize its ii norm, i.e., 

IN|i= EVKIl2 + l|bfcll2- (12) 

fc=0 

According to the preceding discussion, we regularize the 
total variation of the amplitude functions and the £i norm of 
the frequency spectrum subject to the perfect reconstruction 
constraint. This is expressed as the optimization problem: 
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Fig. 1. Example 1: Test signal synthesized as the sum of three sinusoids 
each with a amphtude/phase discontinuity. 



K-1 



mm 

a,b 



^ (lIDafclK + IIDbfelK + Av/|K||^ + ||bfe||^) 

k=0 



K-1 



s. t. x{n) = (^a{n, k) c{n, k) + h{n, k) s{n, k)j (PO) 

fe=0 

where A > is a regularization parameter, and where a, b e 
R^""^ with a = (a(n,fc))„ez,,,,feezK. and afc,bfe € 
with a^; = (a(fc, n))„gZjv The c{n,k) and s{n,k) denote the 
cosine and sine terms, c{n,k) = cos(27r/fcn) and s{n,k) = 
sin(27r/fen), where fk = k/K. 

Notice in (PO), that A adjusts the relative weighting between 
the two regularization terms. If K ^ 2N, then as A — )• 0, 
the minimization problem leads to a solution approaching 
||Da;;||i = ||Dbfe||i = for all frequencies k G I^k- In this 
case, a/fc and bfc are non-time-varying. When K = 2N, they 
coincide with the DFT coefficients, specifically, K[ak{n) — 
']bk{n)) is equal to the DFT coefficient X{k) in (2). 



A. Example 1: A Synthetic Signal 

A synthetic signal x{n) of length A'' = 100 is illustrated 
in Fig. 1. The signal is formed by adding three sinusoidal 
components, Xi{n), i = 1,2,3, with normalized frequencies 
0.05,0.1025 and 0.1625 cycles/sample, respectively. However, 
each of the three sinusoidal components possess a phase 
discontinuity ('phase jump') at n = 50, 38, and 65, respectively. 
The first component, Xi{n), has a discontinuity in both its 
amplitude and phase. 
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Fig. 2. Example 1: Signal decomposition using sparse frequency analysis 
(SEA). The frequency spectrum is sparse and the recovered sinusoidal 
components accurately retain the amplitude/phase discontinuities. 



Here, we set K — 100, so the uniformly-spaced frequencies 
fk, k e Zx, from to 0.5, are separated by 0.005 cy- 
cles/sample. The frequency grid is similar to that of a 200-point 
DFT (including zero-padded). The component xi lies exactly 
on the frequency grid, i.e. 0.05 = 10x0.005. On the other hand, 
the frequencies of components X2 and xz he exactly halfway 
between frequency grid points, i.e., 0.1025 = 20.5 x 0.005 
and 0.1025 = 32.5 x 0.005. The frequencies of X2 and X3 
are chosen as such so as to test the proposed algorithm under 
frequency mismatch conditions. 

Using the iterative algorithm developed below, we obtain 
a{n, k) and h{n, k) solving the optimization problem (PO). The 
frequency spectrum, Zk, defined by (11), is illustrated in Fig. 2. 
The frequency spectrum is clearly sparse. The component Xi 
is represented by a single line in the frequency spectrum at 
k = 10. Note that, because the components X2 and X3 of the 
synthetic test signal are not aligned with the frequency grid, 
they are each represented by a pair of lines in the frequency 
spectrum, at fc = (20,21) and k = (32,33), respectively. This 
is similar to the leakage phenomena of the DFT, except here, 
the leakage is confined to two adjacent frequency bins instead 
of being spread across many frequency bins. 

To extract a narrow-band signal from the proposed decom- 
position, defined by the arrays (a, b), we simply reconstruct 
the signal using a subset of frequencies, i.e.. 



gs(n) — (^a(n, k) c{n, k) + b{n, k) s{n, /c)^ 
kes 



(13) 




0.2 0.3 
Frequency (normalized) 



0.5 



20 



40 



60 



80 



100 



20 



20 



Intrinsic mode function 1 



40 



60 80 



100 



Intrinsic mode function 2 



40 



60 80 



100 



1 


-1 



y2(n) 



20 



40 



60 



80 



100 



-1 



AAAAAAAAAa-aAAAa 



20 



40 60 
Time (samples) 



80 



100 



Fig. 3. Example 1: Signal components obtained by LTI band-pass filtering. 
The Fourier transform X{f) is not sparse and the filtered components do not 
retain the amphtude/phase discontinuities. The amplitude/phase discontinuities 
are spread across time and frequency. 



where S C I^k is a set of one or several frequency indices. 
With S = {10}, we recover a good approximation of Xi. 
With S = {20,21} and S = {32,33}, we recover good 
approximations of components X2 and 2:3, respectively. These 
approximations are illustrated in Fig. 2. Note, in particular, 
that the recovered components retain the amplitude and phase 
discontinuities present in the original components Xi{n). In 
other words, from the signal x{n), which contains a mixture 
of sinusoidal components each with amplitude/phase discon- 
tinuities, we are able to recover the components with high 
accuracy, including their amplitude/phase discontinuities. 

Let us compare the estimation of components Xi from x 
using the proposed method with what can be obtained by band- 
pass filtering. Band-pass filtering is widely used to analyze 
components of signals, e.g., analysis of event-related potentials 
(ERPs) by filtering EEG signals [18], [31]. By band-pass 
filtering signal x using three band-pass filters, Hi, i = 1,2,3, 
we obtain the three output signals, yi, illustrated in Fig. 3. 
The band-pass filters are applied using forward-backward 
filtering to avoid phase distortion. Clearly, the amphtude/phase 
discontinuities are not well preserved. The point in each output 
signal, where a discontinuity is expected, exhibits an attenuation 
in ampUtude. The amphtude/phase discontinuities have been 
spread across time and frequency. 

The frequency responses of the three filters we have used 
are indicated in Fig. 3, which also shows the Fourier transform 
of the signal x{n). Unlike the frequency spectrum in 
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Fig. 4. Example 1: Signal decomposition with the empirical mode 
decomposition (EMD). The first three intrinsic mode functions (EVlFs). The 
amplitude/phase discontinuities are not preserved. 



Fig. 2, the Fourier transform in Fig. 3 is not sparse. 

We also compare SFA and band-pass filtering with empirical 
mode decomposition (EMD) [30]. EMD is a data-adaptive 
algorithm that decomposes a signal into a set of zero-mean 
oscillatory components called intrinsic mode functions (IMFs) 
and a low frequency residual. For the test signal (Fig. la), the 
first three IMFs are illustrated in Fig. 4. (The EMD calculation 
was performed using the Matlab program emd . m from http: 
//perso. ens- lyon.fr/patrick.flandrin/emd.html [52] .) 

Note that the IMFs do not capture the amplitude/phase 
discontinuities. This is expected, as EMD is based on differ- 
ent asstunptions regarding the behavior of the narrow-band 
components (smooth instantaneous amplitude/phase functions). 

Comparing the sparse frequency analysis (SFA) results in 
Fig. 2 with the band-pass filtering results in Fig. 3 and the 
EMD results in Fig. 4, it is clear that the SFA method is better 
able to extract the narrow-band components while preserving 
abrupt changes in amphtude and phase. 



B. Optimization Algorithm 

In this section we derive an algorithm for solving problem 
(PO). We use the alternating direction method of multiphers 
(ADMM) and the majorization-minimization (MM) method to 
transform the original problem into a sequence of simpler 
optimization problems, which are solved iteratively, until 
convergence to the solution. 

By variable splitting, problem (PO) can be written as: 



mm 

a,b,u,v 



K-1 

E 

fc=0 



I Da 



IDb 



fell! 



K-1 
fe=0 
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K-1 



s. t. x{n) = (^u(n,k) c{n,k) + v{n,k) s{n,k)j (14a) 



N-l 



fe=o 
u = a 

V = b 



(14b) 
(14c) 



where u, v G R^^-^ correspond to matrices a and b. 
The augmented Lagrangian [2] is given by 

io(a,b,u,v, A,^) 



fe=0 



,|Dbfc||i+A^K||^ + |K||^ 



+ M ||ufe - a/c - pfella + /i - bfe - qkWl j (15) 

where /i > is a parameter to be specified, and where the 
equality constraint (14a) still holds. The parameter /j. does 
not influence the minimizer of cost function (PO). Therefore, 
the solution of (15) leads to the solution of (PO). The use of 
ADMM yields the iterative algorithm: 

K-l 

a,b-<- argmin ( ||Dafe||i + ||Dbfe||i 

a,b ^— ' \ 

fc=0 

+ ||uft - afe - pfc||2 + H ||vfe - bfe - qfe||2 ) (16a) 



K-l 



u, V are; mm 



|2 II ||2 
\\^k\\2 + l|Vfc||2 



k=0 

+ n ||ufc - afe - pfe||2 + /u ||vfe - bfe - qfe||2 ) 

K-l 

s. t. x{n) = ^2 (^('^j (^{i^i ^) + '^^(^j ^) *(^) (16b) 



fe=0 



P ^ p - (u - a) 
q ^ q - (v - b) 
Go back to (16a). 



(16c) 
(16d) 
(16e) 



In (16a), the variables a and b are uncoupled. Furthermore, 
each of the K columns afe and bfe of a and b are decoupled. 
Hence, we can write 

afe ^ argmin ||Dafe||i + fi ||ufe - afe - pfe||2 (17a) 



bfe argmin ||Dbfe||i +/u||vfe - bfe - qfe| 



(17b) 



for k e Ijk- Problems (17) are recognized as ./V-point TV 
denoising problems, readily solved by [16]. Hence we write 



afe tvd(ufe - pfe, l/yu) 
bfe <(- tvd(vfe - qfe, I 111) 

for k G Ijk- 

Problem (16b) can be written as: 



(18a) 
(18b) 



K-l 

u, V ^ arg mm 

fe=o 

N-l 



N-l 



. J2 \u{n,k)\^ + \v{n,kW^ 

\ n=0 



nm > 

u,v ^ 

fe=0 

N-l 

+ /Lt \u{n, k) — a{n, k) — p{n, k)\ 



n=0 



K-l 



s. t. x{n) = (u{n,k) c{n,k) + v{n,k) s{n,k)j (19) 



fc=o 



which does not admit an explicit solution. Here we use the 
MM procedure for solving (19), i.e., (16b). First we need a 
majorizer. To that end, note that for each k G Z^, 



N-l 



. ^\u{n,k)\^ + \v{n,kW 

\ n=0 



N-l 



^ :;7w E + Hri,k)\') + Ja« (20) 



9 A (») 

^^ife n=0 



where 



N-l 



. ^ |MW(n,fc)|2 + |t;(')(n,A;)|2. (21) 

\ n=0 

Here i is the iteration index for the MM procedure and the 
right-hand side of (20) is the majorizer. An MM algorithm for 
solving (19) is hence given by: 



K-l 



arg mm 



min > 
u,v ^ 



A 



N-l 

J2l-^{\u{n,k)\' + \v{n,kW 

n=0 ^ 



k=0 

+ fj, \u{n, k) — a{n, k) — p{n, k)f 
+ fi \v{n, k) — b{n, k) — q{n, k)\'^ 

K-l 

s. t. x{n) = (^u{n,k)c{n,k) +v{n,k) s{n,k)^. (22) 
fe=o 



The majorizer is significantly simpler to minimize. With the 
use of the quadratic term to majorize the ^2 -norm, the problem 
becomes separable with respect to n. That is, (22) constitutes 
A'^ independent least-square optimization problems. Moreover, 
each least-square problem is relatively simple and admits an 
exphcit solution. Omitting straightforward details, u^'+^' and 
^(i+i) solving (22) are given explicitly by: 



(n, k) = Vk a{n, k) c(n, k) 

+ 2/i(a(n, k) + p{n, k)) 
v^'+'^\n, k) = Vk \a{n, k) s{n, k) 



where 



+ 2n{h{n, k)+qin,k)) 
"-^ _i K-l 



(23a) 



(23b) 



K-l 



fc=0 



fe=0 



7(n, A;) = Vfe c{n,k){a{n,k) + p{n,k)) 



n=0 
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Fig. 5. Example 2: EEG signal, sparse frequency spectrum obtained using 
sparse frequency analysis (SFA), band-pass components reconstructed from 
SPA decomposition. 



+ s(n, k) (b{n, k) + q{n, k)) 



and 



Vk 



A 



2Mfe ' + A 



(24) 



(25) 



for n G Zjv, k G 'Lk- Equations (21) and (23) constitute an 
MM algorithm for computing the solution to (16b). Numerous 
MM iterations are required in order to obtain an accurate 
solution to (16b). However, due to the nesting of the MM 
algorithm within the ADMM algorithm, it is not necessary to 
perform many MM iterations for each ADMM iteration. 

C. Example 2: EEG Signal 

An electroencephalogram (EEG) is the measurement of 
electrical activity of the brain. Several frequency bands have 
been recognized as having physiological significance, for 
example: 4 ^ / ^ 8 Hz (theta rhythms), 8 ^ / < 13 Hz 
(alpha rhythms), 13 < / < 30 Hz (beta rhythms). These 
rhythms are usually obtained by band-pass filtering EEG data 
[50]. 

A one-second EEG signal from [50], with sampling rate 
fs = 100 Hz, is illustrated in Fig. 5. In this example, we aim 



(d) y3(n) 
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Fig. 6. Example 2: Fourier transform of EE signal and band-pass components 
obtained using LTI band-pass filtering. 



to estimate the three noted rhythms via the proposed sparse 
frequency analysis method, and compare the result with that 
of band-pass filtering. 

The proposed method yields the sparse frequency spectrum 
illustrated in Fig. 5. In order to implement (non-linear) band- 
pass filtering using the proposed method, we simply reconstruct 
the signal by weighting each frequency component by the 
frequency response of a specified band-pass filter, 

9H{n) = ^ \H{fk)\(a{n,k)c{n,k) + b{n,k)s{n,k)y 

This generalization of (13) incorporates the frequency response 
of the band-pass filter, H. 

Three band-pass filters Hi, i = 1,2,3, are designed 
according to the theta, alpha, and beta rhythms. Their frequency 
responses are illustrated in Fig. 5, overlaid on the sparse 
spectrum obtained by the proposed technique. The three 
reconstructed components g^. are illustrated in Fig. 5. It 
can be seen that each component has relatively piecewise- 
constant amplitudes and phases. For example, the component 
gui in Fig. 5(c) exhibits an amplitude/phase discontinuity 
at about t = 0.5 seconds. The component in Fig. 5(d) 
exhibits an amplitude/phase discontinuity at about t = 0.35 
seconds and shortly after i = 0.6 seconds. The (instantaneous) 
amplitude/phase functions are otherwise relatively constant. 

The signals obtained by LTI band-pass filtering are shown 
in Fig. 6. The utilized band-pass filters are those that were 
used for the sparse frequency approach. The Fourier transform 
of the EEG signal is shown in Fig. 6(a). 

Comparing the estimated theta, alpha, and beta rhythms 
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obtained by the two methods, shown in Figs. 5 and 6, it 
can be seen that they are quite similar. Hence, the proposed 
method gives a result that is reasonably similar to LTI band- 
pass filtering, as desired. Yet, the proposed approach provides 
a potentially more accurate estimation of abrupt changes in 
amplitude and phase. In this example, the true components are, 
of course, vmknown. However, the sparse frequency analysis 
approach provides an alternative to LTI filtering, useful in 
particular, where it is thought the underlying components are 
sparse in frequency and possess sparse amplitude and phase 
deriviatives. 



IV. Sparse Frequency Approximation 

In appUcations, the available data y{n) is usually somewhat 
noisy and it may be urmecessary or undesirable to enforce the 
perfect reconstruction constraint in (PO). Here we assume the 
data is of the form y = x + w, where w denotes additive 
white Gaussian noise. In this case, a problem formulation more 
suitable than (PO) is the following one, where, as in basis 
pursuit denoising (BPD) [14], an approximate representation 
of the data is sought. 



K-l 



mm 

a,b 



^ (lIDafelli + IIDbfelli + X^\\eik\\l + \\hk\\l) 

k=0 

N-1 r K-l 

-|-Ai y(n)— (^a{n,k)c{n,k)+b{n,k) s{n,k)^ 



n=0 



fe=o 



(PI) 



The parameter Ai > should be chosen according to the noise 
level. In the following, we derive an algorithm for solving 
(PI). 

Applying variable splitting, as above, (PI) can be rewritten: 



K-l 



mm 

a,b 



"'^ fe=0 

N-1 r K-l 

Ai y{n) — (u{n, k) c{n, k) + v{n, k) s{n, k)j 




s. t. 



As above, we apply ADMM, to obtain the iterative algorithm: 

K-l 

)afc||i + ||DbJ|i 



a, b arg min 

a.b 



fc=0 

+ n \\uk - afe - pkWl + IJ, ||vfc - bfc - qfe||2 ) 

K-l 



u, V <— arg mm 



lUfello + llVfello 



(26a) 
(26b) 



fc=0 



+ n \\uk - aft - TpkWl + M l|vfc - bfc - Qfella ) 

N-lr K-l 

+ Ai y{n) — k)c{n, k) + v{n, k)s{n, kfj 



n=0 



fe=0 



p ^ p - (u - a) 
q <(- q - (v - b) 

Go to (26a). 



(26c) 
(26d) 

(26e) 



Note that step (26a) is exactly the same as (16a), the solution 
of which is given by (18), i.e., TV denoising. To solve (26b) 
for u and v, we use MM with the majorizer given in (20). 
With this majorizor, an MM algorithm to solve (26b) is given 
by the following iteration, where i is the MM iteration index. 



Af-l 



<— arg mm 



nin > 
u,v ^ 



K-l 



^--^{\uin,k)\' + \v{n,k)\' 



n=0 I fc=0 

K-l 



2A 



+ Ai ( y{n) — {u{n, k) c(n, k) + v{n, k) s{n, k)) 
^ fe=o 

K-l 

+ II \u{n, k) — a{n, k) — p{n, A;)^ 



fc=0 
K-l 



+ fj, \v{n, k) — b{n, k) — q{n, 



fe=o 



(27) 



where AJ^ is given by (21). As in (22), problem (27) decouples 
with respect to n and the solution be found 

in expUcit form: 



(n, k) = Vk 



v^'+^\n,k) = Vk 
where 



/?(n, k) c(n, k) + 2^(a(n, k) + p{n, fc)) 
(3{n, k) s{n, k) + 2/x(6(n, k) + q{n, k)) 



^ K-l _^ K-l 



7(n, k) = Vk c{n, k) (a(n, k) + p{n, k)) 

+ s{n, k){b{n, k) + q{n, k)) 
where Vk is given by (25). 



(28) 



A. Example 3: A Noisy Multiband Signal 

This example illustrates the estimation of a sinusoid with a 
phase discontinuity when the observed data includes numerous 
other components, including additive white Gaussian noise. 
The example also illustrates the estimation of the instantaneous 
phase of the estimated component. The resulting instantaneous 
phase function is compared with that obtained using a band- 
pass filter and the Hilbert transform. The Hilbert transform is 
widely used in appUcations such as EEG analsys [38], [46], 
[58], [60]. 

The one-second signal x{t), with sampling rate of Fg = 100 
Hz {N = 100 samples), illustrated in Fig. 7c, is synthesized 
as follows: adding the 9.5 Hz sinusoid r{t) (with phase 
discontinuity at t = 0.4 seconds) in Fig. 7a, and the sinusoids 
and white Gaussian noise illustrated in Fig. 7b. Note that the 
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(a) Sinusoid r(t) witti phase discontinuity 



0.2 0.4 0.6 0.8 

(b) Additional sinsoids and noise 




0.4 0.6 
Time (sec) 

Fig. 7. Example 3: The signal x{t) in (c) is synthesized as the sum of (a) 
9.5 Hz sinusoid with a phase discontinuity and (b) additional sinusoids and 
white Gaussian noise. 



instantaneous frequency of r{t) has an impulse at t = 0.4 due 
to the phase discontinuity. 

The sparse frequency analysis approach (PI) was solved 
with the iterative algorithm described above. We used K = 50 
uniformly spaced frequencies from to 50 Hz, we obtain 
a, b G M^^^, the time- varying cosine and sine amplitudes, 
with discrete frequencies fk = k Hz. The sparse frequency 
spectrum is illustrated in Fig. 8(a). Our aim is to recover r{t), 
but (by design) its frequency of 9.5 Hz is halfway between the 
discrete frequencies /g and /lo, which are clearly visible in 
the sparse spectrum. Hence, an estimate of r{t) is obtained via 
(13) using the index set S = {9, 10}, i.e., f{t) = 5'{9,io}(t), 

9s {t) = ag{t) COSWgt + bg{t) sincogt 

+ aio{t) cos wioi + bio (t) sin wio* • (29) 

The time-varying amplitudes, a,k{t) and for fc = 9, 10, are 
illustrated in Fig. 8 (dashed hues). Within the same plots, the 
functions ak{t) cos{uikt) and bk{t) sm{oJkt) with Wfc = 27r/fc 
for k — 9, 10 are also shown (solid lines). The piecewise- 
constant property of ak{t) and bk{t) is clearly visible. The 
signal gs{t) is also illustrated in the figvffe. Note that gsit) 
has a center frequency of 9.5 Hz, while the sinusoids from 
which it is composed have frequencies 9 and 10 Hz. 

To compute and analyze the instantaneous phase of the signal 
gs{t), it is convenient to express gs{t) in terms of a single 
frequency (9.5 Hz) instead of two distinct frequencies (9 and 
10 Hz). Therefore, by trigonometric identities, we write 
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Fig. 8. Example 3: Signal decomposition using sparse frequency analysis 
(SEA). The discontinuity in the instantaneous phase of the 9.5 Hz sinusoid is 
accurately recovered. 



gs{t) = am{t) cos(w„i) + b^it) sin(a;„i) 

+ am+i (i) cos{u)m+it) + bm+1 {t) sm{u)m+it) (30) 



as 



where 



gs {t) = um (t) cos(wMi) + bM (t) sin(wMi) (3 1) 



auit) = {am{t) + ajn+i{t)) cos(Aut) 

- {bm{t) - bm+i{t)) sin(Awi) (32) 
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6m (i) = {bm{t) + bm+i{t)) cos(Awt) 

- {am{t) - am+i{t)) sin(Awf) (33) 

and 

= + a;m+i)/2, Aw = - w„)/2. 

Here lum is the frequency midway between ojm and oJm+i- 
Equation (31) expresses gs{t) in terms of a single center 
frequency, uim, instead of two distinct frequencies as in 
(30). Note that aM{t) and bM{t) are readily obtained from 
the time varying amplitudes ak{t) and bk{t). The functions 
aM(i) cos(cjMi) and bM{t) sin{(jjMt) are illustrated in Fig. 8, 
where aM{t) and bM(t) are shown as dashed lines. Note that 
these ampUtude functions are piecewise smooth (not piecewise 
constant), due to the cos{Aujt) and sin(Aa;t) terms. 

To obtain an instantaneous phase function from (31) it 
is furthermore convenient to express gs{t) in terms of 
M{t) exp(j wm*)- To that end, we write gs{t) as 

+ ^bM (t) ei""* -^bM (i) e-J""* (34) 



or 



9s{t) 



+ 



which we write as 



(35) 



(36) 



where g+{t) is the 'positive frequency' component, defined as 



g+{t) := auit) + j buit). 



(37) 



According to the model assumptions, g+{t) is expected to be 
piecewise smooth with the exception of sparse discontinuities. 
We can use (37) to define the instantaneous phase function 



7M 



(t) 



tan 



_fbM{t)\_ 
\aM{t))' 



(38) 



The function 6m (t) represents the deviation of gs (t) around 
its center frequency, wm- It is useful to use the four-quadrant 
arctangent function for (38), i.e., 'atan2'. For the current 
example, the instantaneous phase 9m (t) for the 9.5 Hz signal, 
gs{t), is illustrated in Fig. 8. It clearly shows a discontinuity 
att = 0.4 of about 180 degrees. 

A finer frequency discretization would also be effective 
here to reduce the issue of the narrow-band signal component 
(9.5 Hz) falling between discrete frequencies (9 and 10 Hz). 
However, an aim of this example is to demonstrate how this 
case is effectively handled when using sparse frequency analysis 
(SFA). 

The estimate of the 9.5 Hz signal r{t) obtained by band- 
pass filtering is illustrated in Fig. 9a (sohd line). The utilized 
band-pass filter was designed to pass frequencies 8-12 Hz 
(i.e., alpha rhythms). The Hilbert transform is also shown 
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(a) Analytic signal via BPF/Hllbert filtering 
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Fig. 9. Example 3: The estimate of the 9.5 Hz component using LTI band- 
pass filtering. The instantaneous phase, computed using the Hilbert transform, 
exhibits a gradual phase shift. 



(dashed line). The instantaneous phase of the analytic signal 
(around 9.5 Hz) is illustrated in Fig. 9b. It can be seen that 
the instantaneous phase undergoes a 180 degree shift around 
t = 0.4 seconds, but the transition is not sharp. Given a real- 
valued signal y{t), a complex-valued analytic signal ya{t) is 
formed by ya{t) = y{t) + j yuit), where ynit) is the Hilbert 
transform of y{t). 

In contrast with LTI filtering (band-pass filter and Hilbert 
transform), the sparse frequency analysis (SFA) method yields 
an instantaneous phase function that accurately captures the 
step discontinuity at i = 0.4. Hence, unlike LTI filtering, the 
sparse frequency analysis (SFA) approach makes possible the 
high resolution estimation of phase discontinuities of a narrow- 
band signal buried within a noisy multi-band signal. This is 
true even when the center frequency of the narrow-band signal 
falls between the discrete frequencies fk- 

V. Measuring Phase Synchrony 

The study of synchronization of various signals is of interest 
in biology, neuroscience, and in the study of the dynamics 
of physical systems [33], [47]. Phase synchrony is a widely 
utilized form of synchrony, which is thought to play a role in 
the integration of functional areas of the brain, in associative 
memory, and motor planning, etc. [33], [57], [59]. Phase 
synchrony is often quantified by the phase locking value (PLV). 
Various methods for quantifying, estimating, extending, and 
using phase synchrony have been developed [5], [6], [28]. 

The phase synchrony between two signals is meaningfully 
estimated when each signal is approximately narrow-band. 
Therefore, methods to measvffe phase synchrony generally 
utilize band-pass filters designed to capture the frequency band 
of interest. Generally, the Hilbert transform is then used to 
compute a complex analytic signal from which the time-varying 
phase is obtained. In this example, we illustrate the use of SFA 
for estimating the instantaneous phase difference between two 
channels of a multichannel EEG. 

A. Example 4: EEG instantaneous phase difference 

Two channels of an EEG, with sampling rate fs = 200 Hz, 
are shown in Fig. 10. Band-pass signals in the 11-15 Hz band. 
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Fig. 10. Example 4: Two channels of a multichannel EEG signal. 
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Fig. 11. Example 4: Band-pass (11-15 Hz) signals estimated using LTI band- 
pass filtering. Instantaneous phase functions obtained using Hilbert transform. 
Channel 1 (left) and channel 2 (right). Bottom: instantaneous phase difference. 

obtained via band-pass filtering are illustrated in Fig. 11. The 

instantaneous phase functions (around 13 Hz) are computed 
using the Hilbert transform, and the phase difference is shown. 
The computation of the phase locking value (PLV) and other 
phase-synchronization indices are based on the difference of 
the instantaneous phase functions. 

The sparse frequency analysis (SFA) technique provides an 
alternative approach to obtain the instantaneous phase of each 
channel of the EEG and the phase difference function. Phase 
synchronization indices can be subsequently calculated, as 
they are when the instantaneous phase difference is computed 
via LTI filtering. Here we use problem formulation (PO) with 
K = 100, with the frequencies fk, equally spaced from to 99 
Hz, i.e., fk = k Hz, ^ k ^ K — 1. With this frequency grid, 
the 11-15 Hz band corresponds to five frequency components, 
i.e. k e S = {11, . . . , 15}. The 11-15 Hz band can then be 
identified as gs{t) in (13). 



In order to compute the instantaneous phase of gs{t), 
we express gs{t) in terms of sins/cosines with time-varying 
amplitudes and constant frequency as in (31). For this purpose, 
we write gs{t) as 

9s{t) = 5{ii,i5}(i) + fi'{i2,i4}(i) + 5{i3}(i)- (39) 
Using (32) and (33), we can write 

^{is} (i) = Om (*) cos{uMt) + h^^ {t) sm{ujMt) (40) 

5{i2,i4} (0 = (*) cos((x;Mi) + b^M (0 sin(wMf) (41) 

5(11,15} {t) = Om (*) cos(a;Mi) + &m {t) sin(u;MO (42) 

where ujm = 27r/i3 = 267r, with /13 being the middle of 
the five frequencies in S. The functions a^^(t), are 
obtained by (32) and (33) from the amplitude functions ak{t), 
hk{t) produced by SFA. Therefore, we can write the 11-15 Hz 
band signal, gs{t), as (31) where 

aM (t ) = aSV) + «M ) + OmV) (43) 

hM{t)^b'^^{t) + h^l^{t)+b^^{t). (44) 

Further, the instantaneous phase of gs{t) can be obtained using 
(38). The 11-15 Hz band of each of the two EEG channels so 
obtained via SFA, and their instantaneous phase functions, are 
illustrated in Fig. 12. 

Comparing the phase difference functions obtained by SFA 
and BPF/Hilbert (LTI) filtering, it can be noted that they are 
quite similar but the phase difference estimated by SFA is 
somewhat more stable. Between 0.3 and 0.4 seconds, the 
SFA phase-difference varies less than the LTI phase-difference. 
Moreover, in the interval between 0.15 and 0.2 seconds, the LTI 
phase difference is increasing, while for SFA it is decreasing. 

The true underlying subband signals are unknown; yet, if 
they do possess abrupt amplitude/phase transitions, the SFA 
technique may represent them more accurately. In turn, SFA, 
may provide more precise timing of phase locking/unlocking. 

VI. Conclusion 

This paper describes a sparse frequency analysis (SFA) 
method by which an A^-point discrete- time signal x{n) can be 
expressed as the sum of sinusoids wherein the ampUtude and 
phase of each sinusoid is a time-varying function. In particular, 
the amphtude and phase of each sinusoid is modeled as being 
approximately piecewise constant (i.e., the temporal derivatives 
of the instantaneous amplitude and phase functions are modeled 
as sparse). The signal x{n) is furthermore assumed to have a 
sparse frequency spectrum so as to make the representation well 
posed. The SFA method can be interpreted as a generalization of 
the discrete Fourier transform (DFT), since, highly regularizing 
the temporal variation of the amplitudes of the sinusoidal 
components leads to a solution wherein the amplitudes are 
non-time- varying. 

The SFA method, as described here, is defined by a convex 
optimization problem, wherein the total variation (TV) of the 
sine/cosine amplitude functions, and the frequency spectrum 
are regularized, subject to either a perfect or approximate 
reconstruction constraint. An iterative algorithm is derived 
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Fig. 12. Example 4: Band-pass (11-15 Hz) signals estimated using sparse 
frequency analysis (SFA). Channel 1 (left) and channel 2 (right). Bottom: 
instantaneous phase difference. 



using variable splitting, ADMM, and MM techniques from 
convex optimization. Due to the convexity of the formulated 
optimization problem and the properties of ADMM, the 
algorithm converges reliably to the unique optimal solution. 

Examples showed that the SFA technique can be used to 
perform mildly non-linear band-pass filtering so as to extract 
a narrow-band signal from a wide-band signal, even when the 
narrow-band signal exhibits amplitude/phase jumps. This is 
in contrast to conventional linear time-invariant (LTI) filtering 
which spreads amplitude/phase discontinuities across time and 
frequency. The SFA method is illustrated using both synthetic 
signals and human EEG data. 

Several extensions of the presented SFA method are envi- 
sioned. For example, depending on the data, a non-convex 
formulation that more strongly promotes sparsity may be suit- 
able. Methods such as re-weighted LI or L2 norm minimization 
[10], [26], [63] or greedy algorithms [36] can be used to address 
the non-convex form of SFA. In addition, instead of modeling 
the frequency components as being approximately piecewise 



constant, it will be of interest to model them as being piecewise 
smooth. In this case, the time-varying amplitude functions 
may be regularized using generalized or higher-order total 
variation [8], [29], [32], [53], or a type of wavelet transform 
[9]. Incorporating higher-order TV into the proposed SFA 
framework is an interesting avenue for further study. 
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